NASA Contractor Report 4203 


Hierarchic Plate and Shell Models 
Based on p-Extension 


Barna A. Szabo and Glenn J. Sahrmann 
Washington University 
St. Louis, Missouri 


Prepared for 

Langley Research Center 

under Grant NAG 1-639 


rvj/NSA 

National Aeronautics 
and Space Administration 

Scientific and Technical 
Information Division 


1988 


TABLE OF CONTENTS 


1. Introduction 1 

2. Conventional and hierarchic theories for plates and shells 4 

2.1. Kirchhoff’s theory 4 

2.2. The Reissner-Mindlin theory 5 

2.3. Higher order theories 6 

2.4. Hierarchic theories 6 

3. Hierarchic beam, arch, plate and shell models 8 

3.1. Beam-columns and arches 8 

3.2. Plates and shells 16 

3.3. Laminated arches and shells 17 

4. Hierarchic basis functions for 19 

4.1 Hierarchic basis functions for S™ 19 

4.2 Hierarchic basis functions for $ p - p -« 20 

5. Examples 23 

5.1. Circular arch, constant cross section, r/t » 15 23 

5.2. Circular arch, constant cross section, r/t ss 1000 25 

5.3. Circular arch, variable cross section 28 

5.4. Cylindrical shell, r/t = 100 30 

6. Summary and conclusions 38 

7. References 40 


PRECEDING PAGB BLANK NOT FILMED 

iii 

9 



1. INTRODUCTION. 


Assurance of the reliability and accuracy of computed data is fundamentally 
important in computer aided analysis and design. In this paper we address the 
question of how to ensure the reliability and accuracy of computed data in engi- 
neering computations concerned with analyses of structures comprised of beam, 
arch, plate and shell components and components which require fully three dimen- 
sional representation. We consider only formulations based on the linear theory 
of elasticity but our approach can be generalized to cases that involve geometric 
and/or material nonlinearities. 

Beams and arches are three dimensional bodies characterized by the fact that 
two of the three dimensions are much smaller than the third. Similarly, plates 
and shells are three dimensional bodies characterized by the fact that one of the 
dimensions is much smaller than the other two. The various theories for beams, 
arches, plates and shells recognize and exploit this. These theories are useful 
because the quantities of interest in the analyses of beams, arches, plates and 
shells, such as membrane forces, bending moments and shear forces, are related to 
certain averages of the displacement across the small dimension(s) of these three 
dimensional bodies. This permits reduction of the dimensions in the case of beams 
and arches from three to one and in the case of plates and shells from three to 
two. 

In engineering problems the line of demarkation between problems of three di- 
mensional elasticity and problems which can be modeled with conventional beam, 
arch, plate and shell theories is not sharp. Usually plates and shells are stiffened 
and/or joined with solid bodies. Of greatest engineering interest are the neighbor- 
hoods of shell intersections, cutouts, attachments, etc. where the stress states are 
truly three dimensional and therefore the assumptions of conventional shell theo- 
ries do not hold. If we are to ensure the reliability and accuracy of computed data 
without sacrificing computational efficency then we must be able to model these 
parts of the structure with three dimensional theories while retaining the simplify- 
ing assumptions incorporated in plate and shell theories where those assumptions 
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hold. This raises the fundamentally important question stated by Naghdi in [1] as 
follows: 

“Under what circumstances do the equations of shell theory supply an ap- 
proximate solution to the three-dimensional equations and how ‘ close * is this 
approximate solution to the exact solution?” 

In this paper we present a systematic process which provides means for find- 
ing answers to this question with respect to specific problems in engineering design 
and analysis. Our approach is based on hierarchic sequences of approximation con- 
structed in such a way that the approximate solutions corresponding to a hierarchic 
sequence of models converge to the exact solution of the fully three dimensional 
model. Selection of the discretization parameters and the stopping criterion are 
based on (1) estimation of the relative error in energy norm; (2) equilibrium tests, 
and (3) observation of the convergence of quantities of interest. This approach 
is closely related to p-extension procedures which have been used successfully for 
estimating and controlling errors of discretization. 

Several beam, arch, plate and shell theories exist. These theories have been 
created and justified by two approaches: 

(a) By a priori assumptions concerning the mode of deformation. This approach 
is favored in the engineering literature, see for example [2-4]. 

(b) By power series expansion of the solution of the three dimensional differential 
equations of elasticity so that powers of the thickness parameter are factored. 
There are several possible variants of this approach: The power series expan- 
sion can be applied to the differential equations of elasticity directly (see for 
example [5-8]) or any of the variational formulations of the differential equa- 
tions of elasticity. Power series expansion procedures applied to variational 
formulations lead to theories characterized by the variational formulation. 
Ciairlet amd Destuynder showed, without a priori assumptions based on phys- 
ical arguments, that Kirchhoff’s theory of plates is the first in a sequence of 
plate theories that can be constructed from the Hellinger-Reissner variational 
principle [9]. 

In this paper we will be concerned with formulations based on the principle 
of virtual work or, equivalently, the principle of minimum potential energy. Our 
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focus is not on the development or justification of new theories for beams, arches, 
plates or shells but on aspects of computer implementation of hierarchic sequences 
of finite element spaces suitable for numerical treatment of a large variety of prac- 
tical problems which may concurrently contain thin and thick plates and shells, 
stiffeners, and regions where truly three dimensional representation is required. 
The discretization parameters which characterize the transverse variation of the 
displacement components are not fixed a priori, but taken into consideration in 
the selection of discretization. 

We have selected the principle of virtual work as the basis for our formulation 
because this is the best understood formulation among alternatives and we have 
substantial experience with it. Analogous construction of hierarchic approximation 
spaces is possible for formulations based on other principles. 

Error control procedures require feedback information concerning the accu- 
racy of the solution in terms of the quantities of interest, and means for reducing 
the error when necessary. Our approach makes it possible to select sequences of 
discretization by adaptive or feedback procedures. We will outline and demon- 
strate such procedures by examples. 
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2. CONVENTIONAL AND HIERARCHIC THEORIES FOR 
PLATES AND SHELLS 


We note that there are fundamental differences between the motivations un- 
derlying the development of classical and modern approaches to modeling of plates 
and shells. Development of classical theories was motivated by the recognition 
that the system of partial differential equations of three dimensional elasticity is 
intractable analytically except in severely restricted cases. Reduction of the num- 
ber of dimensions in the case of beams, arches and bars from three to one and 
in the case of plates and shells from three to two permit analytical treatment of 
large classes of problems. In the cases of arches and shells the coordinate systems 
must be appropriately chosen (e.g. cylindrical, spherical, ellipsoidal, etc. systems) 
to allow analytical treatment. Comprehensive surveys of classical theories with 
historical notes and lists of key references are available in [1,10]. 

The motivation of modern development is quite different. The main goal is 
to allow computer implementation so that a very wide range of problems can be 
analyzed by numerical methods efficiently and with guarantee of reliability. The 
range of problems is to include, for example, simple bars as well as laminated shells 
of arbitrary curvature, regions of shell intersections and solid bodies attached to 
shells. 

In the following we briefly review the essential features of the most widely 
used conventional plate and shell theories and outline the hierarchic theory. The 
notation used for representing the components of the displacement vector « is 
given in Fig. 2.1. 

2.1 Kirchhoff’g theory. 

In Kirchhoff’s theory of plates, formulated in 1850 [ll], there are three dis- 
placement fields. The functions u„ 0 (*. y), u v0 (x, y) in the following equations repre- 
sent the components of the in-plane displacement vector in the x and y directions, 
respectively, and the function u t0 {x,y) represents the transverse displacement vec- 
tor component. 

«* =u l0 (z, y) - x (2.1a) 
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Fig. 2.1. Components of the displacement vector. Notation. 

=tt v o(*,y) ~ z ~jr~ ( 2 - 16 ) 

dy 

u,=«,o (*,»)• ( 2 - lc ) 

In its generalization to shells, known the Kirchhoff-Love theory [12], u,, u v , «, are 
understood to mean the curvilinear (contravariant) displacement components. 

From the assumed mode of deformation: e, d = f du ,/dz = o. The assumption 
that <r a = 0 implies that the plate is orthotropic with Poisson’s ratios v x , = v v , = 0. 
The shear strains 7 *, and 7 V , are zero. E.g.: 

dug du g du g Q dugQ 

' 7l * ~ ~d7 + ~a^~ dx + dx~ 

In applications of the principle of virtual work the transverse displacement «,o 
and its first derivatives must be continuous. This is a major disadvantage of this 
formulation because enforcement of slope continuity is difficult in the general case. 

2.2. The Reissner-Mindlin theory. 

In the Reissner-Mindlin theory, formulated in the 1940’s and early 1950’s [13- 
15], two new fields «.i, u yl , are introduced. These fields represent the rotation of 
the cross-sectional planes to which (respectively) the x-axis and y-axis is normal: 

u* =u s q[x, y) - zu x i{x, y) (2.2a) 

u v =«vo(*, y) - z«»i(z> y) (2.26) 

u, =u M o(x, y) (2.2c) 
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Again, c, = 0 and a, — 0 is assumed. The shear strains 7 „ and 7 y , are independent 
of z. For example: 

/ . . du„o 

'ix* = y ) + -fo- 

In applications of the principle of virtual work the shear strain energy is adjusted 
by a shear factor. In its generalization to shells «*, u y , u, are understood to mean 
the curvilinear (contravariant) displacement components [16]. 

2.S. Higher order theories. 

In higher order theories the displacement fields are typically approximated by 
expressions of the form: 

n 

«» = /<(*) u x< ( x, y) (2.3a) 

i= 0 

n 

= ^2 U (*) «v< (*. y) (2.S6) 

i=0 

m 

u t =^rfi{z)u Mi (x,y). (2.3c) 

«=o 

Usually fi(z) = z i . The various theories differ in the chice of n and m and the 
constitutive law. For example, in the plate theory proposed by Lo, Christensen 
and Wu [17,18] n = 8, m = 2 and the constitutive law is the stress-strain law of 
isotropic elasticity. Discussion of other higher order theories is available in [17]. 

2.4. Hierarchic theories. 

A hierarchic theory is essentially a system of progressively higher order the- 
ories based on the same generalized formulation and the same constitutive law. 
Each theory within a hierarchic system is embedded in all higher order theories in 
that system and the approximate solutions corresponding to progressively higher 
order theories converge to the exact solution of the generalized formulation of the 
fully three dimensional problem. Various hierarchic theories can be constructed 
by choosing alternative generalized formulations, e.g. the Hellinger-Reissner prin- 
ciple, the principle of virtual work, etc. 

In this paper we present a hierarchic system of theories based on the principle 
of virtual work for homogeneous plates. The general form of approximation is the 
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same as in (2.3a,b,c). We propose hierarchic basis functions for the displacement 
components y), u vi (x,y), u Mi [x,y) so that shell elements can be readily joined 
with three dimensional finite elements. The formulation can be extended to apply 
to laminated plates and shells as well. This is briefly discussed in Section 3.3. 
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3. HIERARCHIC BEAM, ARCH, PLATE AND SHELL MODELS. 

In this section the procedures for the computation of stiffness matrices and 
load vectors for hierarchic plate and shell models are described. We begin with the 
simplest representative case for this class of problems, the case of beam-columns 
and arches. 

S.l. Beam-columns and arches. 

Let us consider the Jfeth element of an arch, shown in Fig. 3.1. The element is 
mapped from the standard quadrilateral element, also shown in Fig. 3.1, by some 
smooth mapping functions: 

*= * (fc) U>»?)» y = »7). (£.>?)€ n.t (3.1a) 

where fi,« is the standard element and the superscript (Jfc) refers to the Jfcth finite 
element, fi*. The inverse mapping is: 


£ = £ (fc) (z,y), ft = t] w (x,y), (x, y) € ft*. (3.16) 

In the following we omit the superscript (Jfc) with the understanding that the dis- 
cussion refers to the Jfcth element. 

In the case of beam-columns and arches one dimension is generally much larger 
than the other two. We will assume that the lines corresponding to constant t) 
values (-1 < £ < l) are in the “long” direction of the curved beam or arch. 

Let us now consider an arbitrary point P, as shown in Fig. 3.1. Point P is 
located at the intersection of two coordinate lines, one corresponding to constant 
£, the other corresponding to constant tj values. We denote the unit vector which 
is tangent to the constant i\ line by e ( and the unit vector which is tangent to 
the constant £ line by e n . Given the mapping (3.1) the normalized covariant basis 
vectors are: 
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Fig. 3.1. Typical arch element. Notation. 



Note that is generally not orthogonal to e n . 


(3.2 b) 


We denote the displacement vector components in the direction of (resp. 
e v ) by u ( (resp. Similarly, we denote the displacement vector components in 
the direction of x (resp. y) by «, (resp. u„). The components u*, u v are related to 
u ( , u„ by the following relationships: 





With this notation (3.3a,b) can be written as: 


{«}(*,») = MMu.t?) ( 3 - 5 ) 

where the definition of the 2x2 matrix \R] is obvious from (3.3a,b). 

We now define space $ p -« to be the span of the following (p+l)($ + l) monomials 
on the standard quadrilateral element: 

i, t, .... e, 

v, £»?> £ a >7> •••, Cn, 


ri q , W, eV, .... 

For example, the spanning set for $ *- 2 is the set of monomials inside the dotted 

lines in Fig. 3.2. 
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Fig. 3.2. Space S* ,a is spanned by the monomial terms inside the dotted lines. 
We denote the dimension of S p,t by n, that is: 

n di^^$ p ' , = (p + l)(g + 1). (3.6) 

We define n basis functions for 5 p -« and denote these basis functions by 
i = 1,2,..., n. The definition of basis functions is given in Section 4. Thus any 
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function u, defined on the standard element fi,t, which lies in S p,q can be written 
in the form: 

u = f^a i N i (tn) (3.7) 

«=i 

where o< are arbitrary real numbers. We remark that in the finite element method 
o< are not completely arbitrary because interelement continuity constraints and 
constraints that represent kinematic (principal) boundary conditions must be en- 
forced. Although the basis functions are defined on the standard domain, and are 
given in terms of the standard variables $ and tj, we can (at least in principle) 
substitute (3.1b) for ( and tj and view the basis functions as functions defined on 

n*. 

We will investigate two approaches: In the first approach the approximations 
to the curvilinear components of the displacement vector «*, u v lie in S™. In 
the second approach the approximations to the Cartesian components of the dis- 
placement vector Vy lie in S p '». The first is the classical approach. Its main 
advantage is that for certain types of mapping analytical solutions can be obtained. 
Also, in this approach different degrees of approximation can be used for the dis- 
placement components. For example, in the Reissner-Mindlin theory e 5 P>1 and 
u, e S p>0 . Note that index q identifies a particular beam or arch theory within the 
hierarchic system of theories whereas index p identifies a particular discretization 
within a hierarchic sequence of discretizations based on p-extension. The second 
approach is better suited for computer implementation because enforcement of 
the appropriate continuity conditions between elements mapped by different map- 
ping functions, a condition which frequently occurs in engineering applications, 
is simpler. This permits the use of a great variety of mapping functions. Also, 
computation of the stiffness matrices and load vectors as well as post-solution 
procedures are somewhat simpler. On the other hand, when Cartesian systems 
are used then the displacement components are not oriented in the long and short 
dimensions in general. Therefore the same degree of approximation has to be used 
for each displacement component. While this increases the computational work 
somewhat, an important advantage is gained: all monomial terms which contain 
t), rf i f are retained in the approximation. 
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When we approximate the curvilinear components we have: 


( «i 1 


f u< 1 

'tfi 

n 2 . 

.. N n 

0 

0 . 

o' 


“1 

^2 

{ J = 







< 

. 

l«n j 

0 

0 . 

.. 0 

Ni 

n 2 . 

• • N m _ 


: 


' «n+m ' 


(3.8a) 


or, in short hand: 


{•}(«.,) = WW 


(3.86) 


where ax, a 2 ..., a n+m are coefficients to be determined from the finite element 
solution. When we refer to a specific element, for example the Ath element, we 

write: (t = 1,2 n+ m) or, simply, (a( fc )}. We denote individual columns of 

[iV] by {Ni}. Thus we can also write (3.8a,b) in the following form: 


n+m 

{"}(«,»?) = X) ( 38c ) 

«=1 

When we approximate the Cartesian components then the expressions are 
similar but we must have n = m: 


W(.,v) = X>w>- ( s - 9 ) 

*=i 

We remark that if we accept the restriction n = m then the two approaches can 
be interpreted as the choice of the space S p,t : When the curvilinear components 
are approximated then {«}(*,„) € S p, » where S p> 9 is spanned by [i?]{iVi}, i = 1,2, ...2n. 
When the Cartesian components are approximated then {«}(*,») e S p>t where 5 Pi? 
is spanned by {Ni}, i = l,2,...2n. We now describe the computation of elemental 
stiffness matrices and load vectors for both cases. 

3.1.1. Element stiffness matrices. 

The strain energy of the Jfcth element is of the form: 

Uk = j JJ ( <T * c *+ <T v c v"t" r *v'f*y) t dx dy = — JJ ([■D]{«}(*,v)) [■^][-^]{ u }(*,y) t dx dy (3.10) 

where <r x , c v , r sv are the stress components, e x , e„, ~t xy axe the strain components, 
t = t(x,y) is the thickness; [D] is a differential operator matrix. In the case of plane 
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elasticity: 


[D] = 


't, • 

* h 

a_ a_ 

L dy dx J 


( 3 . 11 ) 


[E] is the material stiffness matrix. For example, in the case of plane stress: 


\E} = 


E 


(1 - *) 


n v 

V 1 
0 0 


0 1 
0 

1 - V 
2 J 


( 3 . 12 ) 


where E and v are, respectively, the modulus of elasticity and Poisson’s ratio. 

We now transform the integral (3.10) so that the integration is on the standard 
element fi, t . We denote the Jacobian matrix of the transformation (3.1a) by [J] 
and its inverse by [J*]: 


[j] 


‘Jll 

J12 

def 

“ dx 
dt 

dy ‘ 

dt 

J21 

J22 m 


dx 

dy_ 




. dr, 

dr, . 


[J ] -1 d = [J*] d = 


'ii 


L^ai 


'13 




( 3 . 18 ) 


and we denote the determinant of the Jacobian matrix by \J\. We define: 


[D*] d = 


Jll dt + Jl2 dr, 


j* A + J* ± 

J2l dC Jn dr, 


J* — + J* — J* — + J* 

7T7 » J 22 ^ J ll~ ' J 12 


di 


dt) 


d* 


d 

dr, 


( 3 . 14 ) 


Clearly, [D*\ is [ D ] written in terms of the variables ( and t,. When we approximate 
the Cartesian components of the displacement vector then we can write (3.10) in 
the following form: 

Uk = lJL (3.15a) 
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The subscript * is a reminder that [ D*\ and |J| are computed from the mapping 
functions of the Jfcth element. When we approximate the curvilinear components 
of the displacement vector then from (3.5) and (3.15a) we have: 

U k = \JJ nt ([2?*][«]{u}«,,)) T [£][£>*][*]{«}«,„) t \J\ dtdr,. (3.156) 

The stiffness matrix for the jfcth element is defined so that: 

Uk = i{a (fc) } T [* (k) ]{aW} (3.16) 

therefore, when we approximate the Cartesian displacement components then the 
elements of are: 

= ll at (m{N i }) T [E][D*){N i }t\J\d^d V . (S.17o) 

When we approximate the curvilinear displacement components then the elements 
of [/£"(*)] are: 

*!?> = ff ([z?*][i2]{Ar < }) r [£:](i?*][/z]{^.}t|j|rf€ < i,. (s.m) 

J J n„ 

These terms are usually computed by numerical quadrature. However when the 
mapping is linear then the integrations can be performed in closed form and the 
speed of computation greatly increased. See, for example, [19,20]. 

3.1.2. Element load vectors: traction loading. 

Let us assume that the load is known in terms of traction components applied 
in the normal and tangential directions on the upper surface of the beam-column 
or arch, that is on the surface corresponding to »j = 1, see Fig. 3.1 and (3.1a). We 
will consider this case only, the other cases are treated analogously. By definition, 
the potential of external traction loads acting on element fl* is: 

P k = I ( T x u x + T v u„) dS (3.18) 

Jo n* 

where d(l k represents the bounding surfaces of element *; T x , T v are components 
of the traction vector; dS is the differential surface. We will write dS in terms of 
the differential arclength do and the thickness: dS = t do ( see Fig. 3.3). In practical 
problems usually T n and T t rather than T x and T v are given. Referring to Fig. 3.3 
we have: 

{£}-[r <»■«> 
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Fig. 3.3. Notation: Traction vector components. 
From the definition of a in Fig. 3.4 we have: 


dx iy 

sin a = — — ; cos a = — 

da da 


and, from (3.1a): 


dx = 


dxW 


dt 


... j d y w 

iy = -af 


n=i 


d(. 


Therefore we can write: 


Pk - £‘ 


dyW dlW 1 

di dz 

dx 5y(*) 


n-i 


dz dz 

If we work in a Cartesian system we substitute: 

or, if we work in the curvilinear system, we substitute: 

= I*] MW 

into (3.22) and write P k in the following form: 

Ph = W T {rW} 


T n 

T t 


tdZ. 


(3.20) 


(3.21) 


(3.22) 


(3.23a) 


(3.236) 


(3.24a) 
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where {rW} is the load vector of the fcth element. Thus, when the Cartesian 
components of the displacement vector are approximated, the »th element of the 
load vector of the fcth finite element is: 



dxW ] 
dyW 

se J 



( 3 . 246 ) 


3.2. Plates and shells. 

The computation of stiffness matrices and load vectors for plate and shell 
elements is analogous to the procedure described in Section 3.1. Plate and shell 
elements are usually six sided but five sided elements are often useful also. We 
will consider six-sided elements only, thus the standard element is the hexahedron, 
shown in Fig. 3.3. The mapping functions are: 

* = *W(£, i 7 , f), y = y (fc) (£>»?>?)> (£,>»,?)» !• (3-25) 

We will assume the mapping to be so that j = 0 is the middle surface of the shell. 
Thus f = +1 gives the ‘upper’ surface and ? = -1 the ‘lower’ surface of the shell. In 
the following we will omit the superscript (Jfc). 




Fig. 3.4. Typical shell element. Notation. 
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The strain energy of the fcth element is: 



= \ Jff mm«}(*.y.')dxdydz 


where (t«}(* >v ,,) d = {«* «„ «,} T is the displacement vector; [D] is the differential 
operator matrix of the three dimensional strain-displacement relationship and [E] 
is the stress-strain law in three dimensions. The definitions are analogous to those 
in (3.11), (3.12) and the details are available in every textbook on elasticity and 
strength of materials. 

Once again we have the option of approximating either the curvilinear or 
Cartesian components of the displacement vector. The relationship between the 
curvilinear and Cartesian components is analogous to (3.5). Thus the approxima- 
tion is either: 

W(e.n.f) = MW (3.27a) 

where [TV] is a 3 x (2» + m) matrix of basis functions; or: 


{«}(.,*.) = MW (3-276) 

where [N\ is a 3 x 3n matrix of basis functions. The basis functions are defined 
on the standard hexahedral element. The space spanned by the basis functions is 
denoted by $ **«. Basis functions, based on Legendre polynomials, are defined in 
the following section. Equation (3.27a) represents the conventional approach. 

3.3. Laminated arches and shells. 

In this paper we are primarily concerned with homogeneous arches and shells. 
In the case of laminated arches and shells the normal and shear stresses are con- 
tinuous across laminar interfaces but the strains are not. Three approaches are 

possible: 

(l) The arch or shell elements can be “stacked” so that the laminar interfaces 
correspond to interelement boundaries. This is feasible when there axe only 
a few laminae, as in the case of arches and shells of sandwich construction . 


- 17 - 


(2) When there are many laminae and interlaminar stresses are of interest then 

the mapping should be so that the laminar interfaces correspond to constant 
ij lines in two dimensions or constant f surfaces in three dimensions. In this 
case Ni({,ri) in two dimensions (resp. in three dimensions) must be 

piecewise polynomials in the t) (resp. j) direction which are continuous at the 
laminar interfaces, with discontinuous derivatives. The derivative ratios are 
computed so that the condition of continuity of the normal and shear stresses 
is satisfied. 

(3) When there are many laminae and the interlaminar stresses are not of interest 
then the material can be treated as homogeneous, and the material properties 
chosen to represent average properties of the laminae. 
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4. HIERARCHIC BASIS FUNCTIONS FOR S™ AND S™-'. 

The key considerations in the selection of basis functions are to ensure nu- 
merical stability and to make the computational effort required for the generation 
of stiffness matrices and load vectors as small as possible. We wish to perform 
p-extensions, therefore the basis functions should be hierarchic, that is the basis 
functions corresponding to polynomial degree p (p > l) should contain, as subsets, 
the basis functions corresponding to polynomal degrees p-1, p-2, . . 1. Compu- 
tational experience has shown that basis functions based on Legendre polynomials 
have good properties from the point of view of numerical stability. We also wish 
to construct our basis functions in such a way that curved beam elements can 
be readily joined with plane elastic elements and plate and shell elements can be 
readily joined with three dimensional finite elements. In the following we define 
the basis functions used for $ p '•* and S p,p,t in the present investigation. 


4.1. Hierarchic basis functions for S p,t . 

The basis functions for S 1 ' 1 are the usual basis functions for four-noded quadri- 
lateral finite elements: 


Ai=j(l-£)(l-„) A s = ±(l + £)(l + , 7 ) 

N 2 =\(1 + S){1-ti) A 4 = j(1-0(1 + »?) 


( 4 . 1 ) 


where the superscript (°) indicates that these basis functions are associated with 
the vertices which are domains of dimension zero. We now define: 


= 



Pi-iit) & 


i = 2, 3, . . . , p 


(4.2) 


where Pi, P 3 , ..., P p _i are the Legendre polynomials. The basis functions for 5 p l , 
p > 1 are comprised of the four vertex modes (4.1) and the following functions 
associated with the sides corresponding to t] = ±1, hence are usually referred to as 
side modes: 


A 2< _ 8 =i-!fc(0, 




* — 2, S 


(4.3) 
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The basis functions for S p - q , g > 1 are comprised of the basis functions for S p>1 plus 
the following edge modes: 

N 2p+ v- i = N 2p+ 2i-2=^-<l>jM, j = 2,3,...,q (4.4) 

and, when p,q> 2, we have the following basis functions, called internal modes: 

N k = {l -t 9 )(l-ri 2 )Pi{t)P i (ri), i = 0,1,2,. ..,p-2; ; = 0, 1, 2 - 2 (4.5) 

where the superscript 2 indicates that these basis functions are associated with the 
area of the standard quadrilateral element which is a domain of dimension two. 
The subscript k ranges from 1 to (p - 1)(? - 1). 

In summary, we have four nodal basis functions, given by (4.1); 2(p-l)+2(g-l) 
edge modes, given by (4.3) and (4.4); (p - 1)(? - 1) internal modes (p, q > 2), given 
by (4.5), a total of (p+ l)(g + 1) basis functions, which is the dimension of S p - q . 
These basis functions are polynomials and are linearly independent. There are no 
monomial terms in the basis functions which are not in the spanning set for S p > q . 
Therefore these basis functions also span $ p - q . 

4.2. Hierarchic basis functions for £?>?>?. 

4.2.1. The vertex modes: Basis functions for 5 1 * 1 * 1 . 

O 

There are 8 vertex modes, denoted by N i} i = 1,2 8. Once again the 

superscript (°) indicates that these modes are associated with the vertices which 
are of dimension zero: 



(4.6a) 

O 1 

tf*«g(l + 0(1 -•?)(!-?) 

(4.66) 


-«?)(! + f). (4.6c) 

4.2.2. Edge modes. 

We use the superscript 1 to indicate that the edge modes pertain to the edges 
which are of dimension one. 
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(a) In the case of S p,Ptl ( p > 2) there are 8(p - 1) edge modes. These are defined in 
terms of the function fa, see (4.2): 


N 1+sii .2)=\Mm - ikw ) 

( 4 . 7a ) 

•^ 2 + 8 (»- 2 ) = J^*(*?)(1 + £)(1 ~ f ) 

( 4 . 76 ) 


( 4 . 7c ) 


where * = 2,5, . ..,p. 

(b) In the case of S p ' Pt1 (p > 2, g > 2 ) there are 8(p - 1) + 4(g - 1) edge modes. These 
are defined as follows: 


^i+. (< -2)=j^(0(l-l)(l-f) 

( 4 . 8a ) 

•Wa+8(i-2) =^4i(v)(l + f)(! ~ f) 

( 4 . 86 ) 

N S+8(i-2) =J^<(»?)(1 “ OU + f) 

( 4 . 8c ) 

•^9+4(i-2) =J^i(f)(l “ f)( 1 “ •?) 

( 4 . 8d ) 

■^ 10 + 40 - 2 ) + ?)(i - n) 

( 4 . 8c ) 

•^ll+4(i-2) =^i(?)(l + OC 1 + *?) 

( 4 . 8 /) 

•^12+4(i-2) = ^j(f)(l - 0(1 + *?) 

( 4 . 8g ) 


where * = 2, 3 p;/ = 2,3 

4.2.3. Face modes. 

We use the superscript 2 to indicate that the face modes pertain to faces 
which are of dimension two. The hierarchic face modes can be treated in one 
of two ways: Either all face modes are analogous to the internal modes (4.5), in 
which case there are 2(p - l) 3 + 4(p - l)(g - 1 ) face modes, p, g > 2, or only the face 
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modes which are nonvanishing on faces £ = ±1; r\ = ±1 are analogous to (4.5) and 
the face modes which are nonvanishing on faces £ = ±1 are defined by: 


^m = ^(l-e 2 )(l-»7 3 )(l±f)^(0^(»7) *’,y = 0,l p-4; (*' + j) = 0, 1, . . . , p — 4 (4.9) 

and m = !»(*, j) depends on the numbering scheme used. In this case the number 
of face modes (n f ) is: 


4(p — l)(g — 1) for p = 2, S 

4(p-l)(?-l) + (p-2)(p-3) for p > 4. 


(4.10) 


This is analogous to the definition of hierarchic internal modes in the case of plane 
elastic elements [21,22]. 


4.2.4. Internal modes. 

We identify the internal modes by the superscript 3. There are (p-3)(p-2)(g- 
l)/2 internal modes (p > 4; q > 2), defined as follows: 

N m = (l~ f 2 )(l - 0(1 - < 2 )Pi(t)P:(l)PkU) (4.Ha) 

where: 

= 0,1,. ...p — 4; (t + /) = 0, 1, . . . , p — 4; * = 0,1,. ...g- 2 (4.116) 

and m = m(i,j,k) depends on the numbering scheme used. 
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5. EXAMPLES. 


The following examples illustrate p-extension procedures and demonstrate 
that conventional curvilinear bases can be replaced by Cartesian bases for a wide 
range of radius to thickness ratios without significant loss in performance. Ex- 
tension procedures are the only means available for the estimation and control of 
errors of discretization. Here we are interested in ensuring that the stress resul- 
tants (i.e. membrane forces, shear forces and moments) are accurate to within one 
percent relative error. 

We report the percent estimated relative error in energy norm, denoted 
by (t r ) E - These estimates are based on the theoretical estimate of error for p- 
extensions. The constants in the theoretical estimate can be computed once the 
strain energies have been computed from finite element solutions corresponding to 
three polynomial degrees p-2, p-1, p; p > 5. Details are available in [21,23,24]. 
The estimated exact values of the strain energy were computed by the computer 
program PROBE [21]. 

In all example problems the mapping functions were generated by the blending 
function method [21]. Thus circular arcs and cylindrical surfaces were represented 
exactly. 

5.1. Circular arch, constant cross section, r/t«15. 

Our first example is the moderately thick circular arch shown in Fig. 5.1. 
The thickness of the arch, i.e. its dimension perpendicular to the x,y plane, is 
unity. It is loaded by parabolically distributed shearing traction in the plane of 
symmetry so that the shear stresses vanish at the top and bottom surfaces and 
F B y = -1.0. Poisson’s ratio is zero. The goal of computation is to determine the 
stress resultants F AX , Fay, Ma o, F B x and M B q. By definition: 

F x = j (cr* cos a + t xv sin a) dA (5.1a) 

where A is the cross sectional area and a is the angle measured from the positive 
x axis to the outward normal to the cross section. Similarly: 

Fy — J ( r xv cos a + <r v sin a) dA (5.16) 
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and 


M° = J [(r* v cos a + <r v sin a) x — (<r x cos a + r xy sin a) y] dA. 


(5.1c) 


The stress resultants were computed from the finite element solutions directly by 
numerical quadrature using twelve Gaussian quadrature points on each boundary 
segment. Of course, the computed stress resultants are not independent. From 
the equations of statics we have: 

Fax = —Fbx > Fay = —Fby = 1.000, Mao = -Mbo- (5.2) 

The exact solution is not known. To obtain a good approximation to the 
exact solution, we solved this as a problem in two dimensional elasticity. The 
computer program PROBE was used [21]. In PROBE the standard polynomial 
space for quadrilateral elements is spanned by all monomial terms of degree p plus 
the terms $ p tj and £»? p for p > 2. The maximal value of p is 8. Details, including 
the definition of element level basis functions, are available in [21,22]. The results 

for p = 1,2 8 are shown in Table 5.1. We see that at p = 8 the equations of 

equilibrium are satisfied to at least four digits. The estimated exact value of the 
strain energy, determined from a sequence of fully two dimensional finite element 
solutions by PROBE , is 0.146660 x 10 a /E where E is the modulus of elasticity. 


Table 5.1. Circular arch, constant cross section, 3 elements, r/t » 15. 
Solution as a problem of elasticity. 


p 

N 

( e <- )e 

3 

44 

14.9 

4 

67 

6.1 

5 

96 

3.6 

6 

131 

2.5 

7 

172 

1.9 

8 

219 

1.5 


Fax Fay 


2.760 1.950 

1.653 1.141 

1.577 0.996 

1.580 1.000 

1.580 1.000 

1.580 1.000 


M a 

Fbx 

-46.17 

3.366 

-28.17 

-1.732 

-25.82 

-1.577 

-25.86 

-1.580 

-25.86 

-1.580 

-25.86 

-1.580 


Fby Mb 


-1.547 51.58 

-0.984 28.07 

-0.996 25.81 

-0.999 25.86 

-0.999 25.86 

-1.000 25.86 


In this problem the estimated error in energy norm is related to the error in 
the average displacement of the loaded cross section. Because by definition the 
energy norm of the error is the square root of the strain energy of the error, 10 
percent error in energy norm corresponds to 1 percent error in the displacement. 
For q = 1 and q = 2 the error in energy norm is virtually constant as p is increased. 
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Fig. 5.1. Circular arch. 

For 5 = 3 there is a substantial decrease of error. This is due to the fact that 
at 5 = 3 representation of the shear deformation is greatly improved. In general, 
if the estimated relative error in energy norm is not changing as p is increased 
while 5 is held constant, then virtually all of the error in energy norm is caused 
by the fact that q is too low. However, as the results in Table 5.2 indicate, large 
error in energy norm does not necessarily mean that the error in stress resultants 
or other quantities of interest is large. The converse of this statement is also 
true, see for example [24]. At (p,q) = (5,1) the estimated relative error in energy 
norm is 7.9 percent, that is 0.079. The relative error of the displacement of the 
loaded cross section is is 0.079 2 = 0.0062, that is 0.62 percent. This level of precision 
is usually more than adequate in elastostatics. In elastodynamics on the other 
hand, computation of the high frequency displacement modes necessitates greater 
precision and thus higher q values. 

5.2. Circular arch, constant cross section, r/twlOOO. 

Our second example is similar to the first, except the arch is now thin: We 
have changed the r 4 from 14.0 to 14.985, see Fig. 5.1. We examine two cases: In 
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Table 5.2 Circular arch, constant thickness, 3 elements, r/t « 15. 

(a) Curvilinear basis. 


p,g 

N 

(c r )E 

Fax 

Fay 

Mao 

Fbx 

Fby 

Mb o 

5,1 

54 

11.2 

1.668 

0.570 

-18.45 

-0.942 

— 1.370 

16.48 

4,1 

46 

7.9 

1.555 

0.969 

-25.30 

-1.544 

-0.991 

25.33 

5,1 

58 

7.9 

1.580 

1.001 

-25.88 

-1.582 

-0.999 

25.88 

6,1 

70 

7.9 

1.581 

1.000 

-25.87 

-1.581 

- 1.000 

25.87 

7,1 

82 

7.9 

1.581 

1.000 

-25.87 

-1.581 

- 1.000 

25.87 

8,1 

94 

7.9 

1.581 

1.000 

-25.87 

-1.581 

- 1.000 

25.87 

5,2 

51 

11.2 

1.668 

0.570 

-18.45 

-0.942 

-1.370 

16.48 

4,2 

69 

7.9 

1.555 

0.969 

-25.30 

-1.544 

-0.991 

25.53 

5,2 

87 

7.9 

1.581 

1.001 

-25.88 

-1.582 

-0.999 

25.88 

6,2 

105 

7.9 

1.581 

1.000 

-25.87 

-1.581 

- 1.000 

25.87 

7,2 

125 

7.9 

1.581 

1.000 

-25.87 

-1.581 

- 1.000 

25.87 

8,2 

141 

7.8 

1.581 

1.000 

-25.87 

-1.581 

- 1.000 

25.87 

5,5 

68 

8.5 

1.684 

0.556 

-18.39 

-0.937 

-1.396 

16.40 

4,5 

92 

2.0 

1.580 

0.952 

-25.25 

-1.541 

-1.025 

25.28 

5,5 

116 

1.6 

1.605 

0.986 

-25.84 

-1.579 

-1.029 

25.84 

6,5 

140 

1.8 

1.600 

0.988 

-25.84 

-1.579 

-1.023 

25.84 

7,5 

164 

1.1 

1.596 

0.990 

-25.85 

-1.580 

-1.016 

25.85 

8,5 

188 

1.0 

1.591 

0.993 

-25.85 

-1.580 

-1.009 

25.85 


(b) Cartesian basis. 


p,q 

N 

{e r )E 

Fax 

Fay 

Mao 

Fbx 

Fby 

Mb o 

3,1 

34 

15.8 

2.761 

1.950 

-46.18 

-3.867 

-1.547 

51.58 

4,1 

46 

8.0 

1.653 

1.141 

-28.18 

-1.733 

-0.984 

28.08 

5,1 

58 

7.9 

1.576 

1.000 

-25.83 

-1.577 

-0.996 

25.82 

6,1 

70 

7.9 

1.580 

1.000 

-25.87 

-1.581 

- 1.000 

25.87 

7,1 

82 

7.9 

1.581 

1.000 

-25.87 

-1.581 

- 1.000 

25.87 

8,1 

94 

7.9 

1.581 

1.000 

-25.87 

-1.581 

- 1.000 

25.87 

3,2 

51 

15.8 

2.761 

1.950 

-46.18 

-3.367 

-1.547 

51.58 

4,2 

69 

8.0 

1.653 

1.142 

-28.18 

-1.733 

-0.984 

28.08 

5,2 

87 

7.9 

1.576 

1.000 

-25.83 

-1.577 

-0.996 

25.82 

6,2 

105 

7.9 

1.580 

1.000 

-25.87 

-1.581 

- 1.000 

25.87 

7,2 

123 

7.9 

1.581 

1.000 

-25.87 

-1.581 

- 1.000 

25.87 

8,2 

141 

7.8 

1.581 

1.000 

-25.87 

-1.581 

- 1.000 

25.87 

3,3 

68 

14.1 

2.784 

1.946 

-46.29 

-3.376 

-1.577 

51.71 

4,3 

92 

2.3 

1.685 

1.132 

-28.28 

-1.741 

-1.021 

28.20 

5,3 

116 

1.6 

1.603 

0.992 

-25.92 

-1.585 

-1.026 

25.92 

6,3 

140 

1.3 

1.603 

0.993 

-25.94 

-1.586 

-1.023 

25.95 

7,3 

164 

1.1 

1.598 

0.995 

-25.93 

-1.585 

-1.016 

25.93 

8,3 

188 

1.0 

1.593 

0.997 

-25.91 

-1.584 

-1.009 

25.91 
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the first case Poisson’s ratio (i/) is zero, in the second case v = 0.3 and plane strain 
conditions are assumed. 

5.2.1. The case i/ = 0. 

The results of computation are shown in Table 5.3. The errors for p < 5 are 
much larger than in the case of the moderately thick arch discussed in Section 5.1. 
Nevertheless, p-convergence is strong and at p = 7,8 the equilibrium relations are 
satisfied to at least three digits of precision. Once again the data show that choice 
of the basis does not affect the accuracy of the solution significantly. For low p 
values the approximate solutions computed with curvilinear bases are worse than 
the approximate solutions computed with Cartesian bases. For high p values the 
curvilinear bases are slightly better, although both approaches yield good results. 

Table 5.3. Circular arch, r/t « 1000, constant cross section, 3 elements, v = 0 . 

(a) Curvilinear basis. 


P>? N ( e r)£ Fax Fay Mao Fbx Fby Mbo 


3,1 

34 

69.2 

4,1 

46 

9.7 

5,1 

58 

0.1 

6,1 

70 

0.0 

7,1 

82 

0.0 

8,1 

94 

0.0 


- 54.812 - 104.234 

- 13.459 - 26.267 

2.035 1.767 

1.642 1.057 

1.610 1.000 

1.610 1.000 


1763.79 131.438 

440.09 28.010 

- 40.27 - 2.662 

- 28.10 - 1.673 

- 27.11 - 1.610 

- 27.12 - 1.610 


- 5.033 - 1969.76 

- 1.549 - 417.16 

- 0.981 42.90 

- 0.999 28.07 

- 1.000 27.11 

- 1.000 27.12 


(b) Cartesian basis. 


p,q N (cr)js Fax Fay Mao Fbx Fby Mbo 


3,1 

34 

70.8 

4,1 

46 

11.2 

5,1 

58 

0.4 

6,1 

70 

0.0 

7,1 

82 

0.0 

8,1 

94 

0.0 


33.982 56.500 

16.860 30.957 

2.351 2.833 

1.565 0.922 

1.611 0.997 

1.612 0.999 


- 988.73 - 90.010 

- 530.15 - 39.502 

- 56.49 - 3.406 

- 25.77 - 1.514 

- 27.09 - 1.608 

- 27.12 - 1.610 


- 0.664 1350.25 

0.903 594.98 

- 0.761 54.06 

- 0.998 25.68 

- 1.000 27.09 

- 1.000 27.12 


The fully two dimensional solution, using the same three-element mesh, was 
computed by PROBE. At p = 8 the following stress resultants were obtained: 
Fax = 1 - 612 ; F AY = 0 . 999 ; M A o = - 27 . 124 ; F B x = - 1 . 610 ; F BY = - 1000 ; M B0 = - 27 . 124 . 
The computed value of the strain energy at p — 8 [N = 219 ) was 0.381787 x 10 */E, 
where E is the modulus of elasticity. 
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5.2.2. The case v = 0.3, plane strain. 

In the case of thin arches and shells Poisson’s ratio does not affect the stress 
resultants significantly but it does affect the displacement and the strain energy. 
It also affects the discretization error in finite element computations. The ap- 
proach presented herein provides for control of the discretization error. This is 
demonstrated in the following. 

If we use the same three-element mesh as before, and q = 2, then we obtain 
the data shown in Table 5.4. We see that the error in stress resultants is large 
at the fixed support. The reason for this is that the radial displacement induced 
by the nonzero Poisson’s ratio is prevented by the fixed support. This excites the 
singularities, causing an oscillatory behavior of the stresses in the element at the 
boundary. We refer to this as the boundary layer effect. To compensate for the 
boundary layer effect we need to introduce a small element which is approximately 

of the same size as the thickness of the arch or shell. The results obtained for a 
four element mesh, which differs from the mesh shown in Fig. 5.1 in that r< = 14.985 

and a fourth element of arc length 0.015 is introduced at the fixed boundary, is 
shown in Table 5.5. We see that there is a marked improvement in convergence 
and the stress resultants are converging to the same values as in the case of v = 0. 

The fully two dimensional solution, using the same four-element mesh, was 
computed by PROBE. At p = 8 the following stress resultants were obtained: 
Fax = 1610; F A y = 1000; M A0 = -27.125; F B x = -1.611; F B y = -1.000; M B0 = -27.128. 
The computed value of the strain energy at p = 8 (N = 295) was 0.347414 x 10 8 /E 
which is almost exactly 0.381787 x 10 8 [l-u 2 )/E, that is (l-t/ 2 ) times the value of the 
strain energy computed for v = 0. Since the strain energy is proportional to the 
applied force times the displacement in the direction of the force, the displacements 
computed by the two methods are also close. 

This example demonstrates that in the case of thin shells very good results 
can be obtained if we replace E with E/( 1 - v 2 ) but otherwise use v = 0. If we 
do not use this simplification then we must take into consideration the boundary 
layer effects in designing the mesh and the value of q has to be at least 2. 

5.3. Circular arch, variable cross section. 

Our third example is the circular arch of variable cross section, shown in 
Fig. 5.2. The r/t ratio ranges from approximately 4 to approximately 15. The 
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Table 5.4. Circular arch, r/t » 1000 , constant cross section, 3 elements, v = 0 . 3 . 

Cartesian basis. 


p ,7 

N 

(c t )e 

Fax 

Fay 

M ao 

Fbx 

Fby 

Mbo 

3,2 

51 

72.0 

40.540 

70.051 

- 1213.92 

- 118.871 

1.883 

1782.86 

4,2 

69 

14.3 

22.197 

43.754 

- 736.68 

- 60.514 

2.459 

909.99 

5,2 

87 

8.0 

0.145 

0.875 

- 14.88 

- 5.073 

- 0.643 

79.05 

6,2 

105 

6.7 

- 1.491 

- 2.574 

42.22 

- 1.476 

- 0.997 

25.11 

7,2 

123 

5.8 

- 1.525 

- 2.503 

41.60 

- 1.614 

- 1.003 

27.17 

8,2 

141 

5.1 

- 1.588 

- 2.522 

42.35 

- 1.608 

- 0.998 

27.09 


Table 5.5. Circular arch, r/t « 1000 , constant cross section, 4 elements, v = 0 . 3 . 

Cartesian basis. 


P ,9 

N 

(e r )E 

Fax 

Fay 

8,2 

69 

71.7 

- 6.295 

3.927 

4,2 

93 

11.5 

0.145 

2.409 

5,2 

117 

0.8 

0.633 

1.557 

6,2 

141 

0.6 

1.155 

1.264 

7,2 

165 

0.6 

1.418 

l.m 

8,2 

189 

0.6 

1.537 

1.043 

9,2 

213 

0.6 

1.584 

1.015 

10,2 

237 

0.6 

1.602 

1.005 

thickness is 

unity. 

To obtain reference 


Mao 

Fbx 

Fby 

Mbo 

- 6.16 

- 119.302 

1.912 

1789.32 

- 34.54 

- 60.552 

2.473 

910.57 

- 27.08 

- 5.051 

- 0.639 

78.72 

- 27.16 

- 1.480 

- 1.001 

25.17 

- 27.14 

- 1.606 

- 1.001 

27.06 

- 27.13 

- 1.611 

- 1.000 

27.13 

- 27.12 

- 1.611 

- 1.000 

27.13 

- 27.12 

- 1.611 

- 1.000 

27.13 


values, the problem was solved as a 


problem 


of elasticity, by PROBE using 3 finite elements. 


The results are shown in Table 


5.6. The estimated exact value of the strain energy, determined from a sequence 
of fully two dimensional solutions by PROBE, is 0.728440 x 10 */E. 


Table 5.6. Circular arch, variable cross section, 3 elements. 
Solution as a problem of elasticity. 


P 

N 

( e r )e 

Fax 

Fay 

M a 

Fbx 

Fby 

M b 

1 

10 

86.8 

0.734 

1.512 

- 22.71 

- 0.394 

- 1.793 

6.08 

2 

27 

63.7 

2.413 

1.105 

- 32.15 

- 2.320 

- 3.431 

34.22 

3 

44 

22.7 

1.694 

1.164 

- 28.88 

- 3.547 

- 2.036 

53.20 

4 

67 

8.7 

1.706 

0.983 

- 26.79 

- 2.185 

- 0.990 

33.89 

5 

96 

4.5 

1.714 

1.000 

- 27.05 

- 1.715 

- 0.951 

27.09 

6 

131 

3.0 

1.714 

1.000 

- 27.04 

- 1.701 

- 0.996 

26.86 

7 

172 

2.3 

1.714 

1.000 

- 27.05 

- 1.712 

- 1.001 

27.03 

8 

219 

1.8 

1.714 

1.000 

- 27.05 

- 1.712 

- 1.001 

27.03 
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s. 



Fig. 5.2. Circular arch, variable cross section. 

The results computed with the curvilinear and Cartesian bases are shown in 
Tables 5.7a, 5.7b. Although in this case the radius to thickness ratio is quite small, 
remarkably good results can be obtained with the lowest order approximation of 
the transverse displacement (q = l). Compare, for example, the results obtained 
with the fully two-dimensional solution in Table 5.6, p = 6 and the results in Table 
5.7b corresponding to ( p , q) = (6, 1). 

For low q values the Cartesian basis is better, for high p values both approaches 
yield similarly good results. Once again, the error in energy norm is large for q = 1, 2 
with a sharp decrease in the error at q = 3. This is due to better representation of 
the shear deformation terms at q > 3. 

5.4. Cylindrical shell, r/t=100. 

Our fourth example is one of the most widely investigated test problems in 
finite element analysis of shells, often called the Scordelis-Lo problem because 
the first finite element analysis of this problem was performed by Scordelis and 
Lo in 1964 [25]. Other finite element solutions were published by Cowper et 
al. [26] and Forsberg [27] in 1970; Dawe in 1975 [28]; MacNeal and Harder in 
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Table 5.7a. Circular arch, variable thickness, 3 elements. 
Curvilinear basis. 



N 

( e r)s 

Fax 

Fay 

M a 

Fbx 

Fby 

Mg 

3,1 

34 

26.0 

1.559 

1.081 

- 26.69 

- 1.209 

- 1.695 


4,1 

46 

23.2 

1.733 

1.032 

- 27.23 

- 1.606 

— 0.955 

25.45 

5,1 

58 

23.2 

1.710 

1.010 

- 26.81 

- 1.711 

- 0.969 

26.96 

6,1 

70 

23.2 

1.711 

1.016 

- 26.89 

- 1.710 

- 0.999 

26.94 

7,1 

82 

23.2 

1.711 

1.015 

- 26.88 

- 1.706 


26.88 

8,1 

94 

23.2 

1.711 

1.015 

- 26.88 

- 1.706 


26.88 

3,2 

51 

15.3 

1.512 

1.072 

- 26.63 

- 1.195 

- 1.663 

19.26 

4,2 

69 

10.1 

1.729 

1.014 

- 27.31 

- 1.614 

- 0.952 

25.62 

5,2 

87 

10.1 

1.708 

0.992 

- 26.90 

- 1.719 

- 0.972 

27.13 

6,2 

105 

10.0 

1.712 

1.000 

- 27.03 

- 1.717 

- 0.999 

27.10 

7,2 

123 

10.0 

1.713 

1.000 

- 27.03 

- 1.713 

- 1.000 


8,2 

141 

10.0 

1.713 

1.000 

- 27.04 

- 1.713 



3,3 

68 

12.2 

1.539 

1.062 

- 26.68 

- 1.189 

- 1.670 

19.17 

4,3 

92 

2.7 

1.740 

1.010 

- 27.35 

- 1.613 

- 0.989 

25.59 

5,3 

116 

2.1 

1.715 

0.989 

- 26.92 

- 1.719 


27.12 

6,3 

140 

1.8 

1.716 

0.999 

- 27.04 

- 1.717 

- 1.022 


7,3 

164 

1.7 

1.715 

0.999 

- 27.04 

- 1.714 


27.04 

8,3 

188 

1.5 

1.713 

1.000 

- 27.04 

- 1.713 



3,4 

85 

12.2 

1.538 

1.062 

- 26.68 

- 1.189 

- 1.670 

19.17 

4,4 

115 

2.6 

1.739 

1.009 

- 27.33 

- 1.613 

- 0.989 

25.59 

5,4 

145 

1.9 

1.713 

0.988 

- 26.90 

- 1.719 

- 1.004 

27.12 

6,4 

175 

1.6 

1.713 

0.998 

- 27.02 

- 1.717 

- 1.021 

27.09 

7,4 

205 

1.4 

1.713 

0.999 

- 27.02 

- 1.714 

- 1.015 

27.04 

8,4 

235 

1.3 

1.712 

1.000 

- 27.03 

- 1.713 

- 1.008 

27.04 
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Table 5.7b. Circular arch, variable thickness, 3 elements. 

Cartesian basis. 


p,q 

N 

( e r)s 

Fax 

Fay 

M a 

Fbx 

Fby 

Mb 

3,1 

34 

19.7 

1.736 

1.428 

- 32.38 

- 8.274 

- 1.840 

49.36 

4,1 

46 

11.4 

1.682 

0.952 

- 26.25 

- 2.037 

- 0.967 

81.74 

5,1 

58 

11.1 

1.716 

1.000 

- 27.06 

- 1.707 

- 0.966 

26.95 

6,1 

70 

11.1 

1.713 

1.001 

- 27.04 

- 1.705 

- 0.998 

26.92 

7,1 

82 

11.1 

1.713 

1.000 

- 27.03 

- 1.712 

- 1.001 

27.02 

8,1 

94 

11.1 

1.713 

1.000 

- 27.03 

- 1.713 

- 1.000 

27.03 

3,2 

51 

19.1 

1.726 

1.424 

- 32.28 

- 3.276 

- 1.838 

49.40 

4,2 

69 

10.2 

1.675 

0.946 

- 26.14 

- 2.036 

- 0.968 

31.78 

5,2 

87 

9.9 

1.712 

0.998 

- 27.01 

- 1.707 

- 0.966 

26.96 

6,2 

105 

9.9 

1.711 

0.999 

- 27.02 

- 1.705 

- 0.998 

26.92 

7,2 

123 

9.8 

1.712 

1.000 

- 27.03 

- 1.712 

- 1.001 

27.02 

8,2 

141 

9.8 

1.713 

1.000 

- 27.03 

- 1.713 

- 1.000 

27.03 

3,3 

68 

16.8 

1.760 

1.418 

- 82.42 

- 3.297 

- 1.852 

49.69 

4,3 

92 

3.3 

1.689 

0.944 

- 26.22 

- 2.047 

- 1.007 

31.89 

5,3 

116 

1.9 

1.722 

0.997 

- 27.08 

- 1.716 

- 0.999 

27.08 

6,3 

140 

1.6 

1.717 

1.000 

- 27.07 

- 1.713 

- 1.022 

27.03 

7,3 

164 

1.4 

1.715 

1.000 

- 27.06 

- 1.719 

- 1.017 

27.11 

8,3 

188 

1.3 

1.714 

1.000 

- 27.05 

- 1.718 

- 1.010 

27.10 

3,4 

85 

16.8 

1.759 

1.418 

- 32.42 

- 3.297 

- 1.851 

49.69 

4,4 

115 

3.2 

1.688 

0.944 

- 26.21 

- 2.047 

- 1.007 

81.89 

5,4 

145 

1.7 

1.719 

0.997 

- 27.05 

- 1.716 

- 0.998 

27.08 

6,4 

175 

1.4 

1.715 

0.999 

- 27.04 

- 1.713 

- 1.021 

27.08 

7,4 

205 

1.2 

1.713 

0.999 

- 27.04 

- 1.719 

- 1.015 

27.11 

8,4 

235 

1.0 

1.713 

1.000 

- 27.04 

- 1.718 

- 1.008 

27.10 
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1985 [29], and Carpenter et al. in 1986 [30]. Displacement data, computed from 
solutions obtained with various finite element computer programs, are summarized 
in [31]. Both displacement data and stress resultants are given in [32], Selected 
data computed from finite element solutions based on shallow shell formulation 
are tabulated in [26] and data computed from finite element solutions based on 
deep shell formulation are given in [28]. Reference data computed from the exact 
solution of the shallow shell formulation are given in [26]. 

5.4.1. Problem statement. 

The shell is shown in Fig. 5.3. It is loaded by its own weight which is 
approximated by uniformly distributed traction of 90.0 lbf/ft 2 acting on the middle 
surface of the shell in the negative z direction. The cylindrical shell is supported 
by diaphragms at the ends. The diaphragms prevent displacement in the * and 
z directions but allow displacement in the y direction. Poisson’s ratio is zero and 
the modulus of elasticity (2?) is S.O x 10® lbf/in 2 . Because there are two planes of 
symmetry, it is sufficient to discretize only one quarter of the shell. In the following 
we refer only to the domain ABCD, see Fig. 5.3. The goal of computation is to 
compute displacements, membrane forces and moments at selected points to about 
one percent relative error. 

5.4.2. Solution as a problem in three dimensional elasticity. 

The problem was solved as a problem in three dimensional elasticity by means 
of PROBE. Four finite elements mapped by quadratic parametric mapping were 
used. Half of the 90.0 lbf/ft 2 traction acting in the negative z direction was applied 
to the upper surface of the shell, half to the lower surface. The mesh and the 
deformed configuration are shown in Fig. 5.4. The computed values of the strain 
energy, the estimated relative error in energy norm ( t r ) E and the displacement 
components (u,) B , (u v ) B are shown in Table 5.8. 

The displacement component (u,) B computed by various computer codes was 
plotted in [31] against N. The exact analytical solution reported for the shallow 
shell formulation is -3.70331 inches [26] and for the deep shell formulation is -3.53 
inches [31]. Our three dimensional result of (a,) B = -3.613 inches agrees to at 
least three significant digits with the results obtained with STARDYNE ’ s QUAD 8 
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Fig. 5.3. Cylindrical shell problem. (Not to scale). 



Fig. 5.4. Cylindrical shell. Deformed configuration. (Quarter model). 

element using 2016 degrees of freedom and differ by less than one percent from the 
results obtained with ABAQUS ’ S4R element which yielded (u,) B = -3.629 inches 
when 2166 degrees of freedom were used [31]. The strain energy of the exact 
analytical solution of the shallow shell formulation is 0.14707 x 10® in • lbf for one 
quarter of the shell [26]. Of course, the exact solution of the three dimensional 
formulation is not known. Our estimate of the exact value of the strain energy 
corresponding to the three dimensional formulation is 0.145049 x 10 5 in • lbf. This 
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Table 5.8. Computed values of the strain energy, estimated 
relative error e r (u rE ) (percent) and displacement components (u a ) B , (u v ) B . 
Solution as a problem in three dimensional elasticity. Four finite elements. 


p 

N 

U(S FE ) 

(in • lbf) 

e r («F£) 

Mb 
M _ 

Mb 

1 

32 

0.942193 X 10 s 

96.69 

-0.081 

-0.003 

2 

104 

0.605197 x 10 4 

76.34 

-0.937 

-0.360 

S 

176 

0.120684 x 10® 

40.98 

-2.927 

-1.525 

4 

300 

0.144160 X 10 8 

7.83 

-3.606 

-1.902 

5 

476 

0.144844 X 10 s 

3.76 

-3.611 

-1.903 

6 

716 

0.144896 x 10 5 

3.25 

-3.611 

-1.903 

7 

1032 

0.144924 x 10 8 

2.94 

-3.612 

-1.903 

8 

1436 

0.144950 X 10 8 

2.61 

-3.613 

-1.904 


value was found by extrapolation based on detailed analyses. The relative errors 
were computed with reference to this value. 

5.4.3. Solution by hierarchic shell model, Cartesian basis, q = 1. 

In the hierarchic model mapping was by the blending function method. Uni- 
formly distributed tractions were applied on the upper and lower surfaces of the 
shell so that the total force acting on the upper surface was the same as the total 
f force acting on the lower surface. The strain energy, the estimated relative error 

in energy norm e r (uF B ), and the displacement components (u,) B , (« y ) B computed 
from solutions obtained by means of the hierarchic shell model, Cartesian basis, 
q = 1 with one finite element (resp. four finite elements) are given in Table 5.9a 
(resp. Table 5.9b). On comparing the results in Tables 5.9a,b with those in Ta- 
ble 5.8, we see that the simplest hierarchic shell model yields results which are 
of similar quality as the results obtained with the fully three dimensional model 
when v — 0. The strain energy of the fully three dimensional model at p = 8, 
N = i486 is slightly lower than the strain energy of the hierarchic shell model, 
Table 5.9b, (p, q) = 8, 1. This result appears to be inconsistent because the finite 
element meshes are the same and 5 8,1 is a subset of the standard polynomial space 
in PROBE corresponding to p = 8. In reality there is no inconsistency, however: 
The finite element space defined by the shell model is not a subset of the finite 
element space defined by the fully three dimensional model. This is because in the 
fully three dimensional model quadratic parametric mapping was used whereas in 
the hierarchic shell model mapping was by the blending function method so that 


I 


Table 5.9a. Computed values of the strain energy, estimated 
relative error e r (u FE ) (percent) and displacement components (u x )b. 

Solution by hierarchic shell model. Cartesian basis. One finite element. 


P,9 

N 

U(u FE ) 

(in • lbf) 


Mb 

(in) 

Mb 

(in) 

4,1 

64 

0.856464 x 10 4 

63.99 

-1.799 

-8.730 

6,1 

126 

0.144514 X 10 5 

6.07 

—3.623 

-1.912 

8,1 

212 

0.144881 x 10® 

3.40 

-3.610 

-1.903 

10,1 

322 

0.144924 X 10® 

2.94 

-3.612 

-1.904 

12,1 

456 

0.144957 X 10 5 

2.52 

-3.613 

-1.904 


Table 5.9b. Computed values of the strain energy, estimated 
relative error e r (ii FE ) (percent) and displacement components (u,) B , («*)b* 
Solution by hierarchic shell model. Cartesian basis. Four finite elements. 


P>$ 

N 

U(u FE ) 

(in • lbf) 

«r(«F£;) 

Mb 

(in) 

(«*)b 

(in) 

4,1 

224 

0.144130 x 10 5 

7.96 

-3.606 

-1.902 

6,1 

456 

0.144903 X 10 8 

3.17 

-3.611 

-1.903 

8,1 

784 

0.144952 X 10® 

2.58 

-3.613 

-1.904 

10,1 

1208 

0.144993 X 10® 

1.96 

-3.615 

-1.905 

12,1 

1728 

0.145022 X 10® 

1.36 

-3.616 

-1.906 


Table 5.10a. Displacement, moment and membrane force data at points B and C. 
Solution by hierarchic shell model. Cartesian basis. One finite element. 


n n 

N 

Me 

m c 

m c 

(My) C 

me 

m B 

me 

P> H 

(in) 

(lbf /in) 

(lbf/in) 


m_ 

(lbf/in) 

_(lbf)_ 

4,1 

64 

0.0668 

-4088.2 

10778.0 

169.4 

1120 

4700 

-539.9 

6,1 

126 

0.5353 

-15325.0 

-72.4 

61.2 

1951 

6408 

-642.2 

8,1 

212 

0.5422 

-304.9 

-117.8 

104.4 

2088 

6328 

-634.8 

10,1 

322 

0.5417 

-297.7 

-130.8 

101.5 

2060 

6306 

-644.3 

12,1 

456 

0.5412 

-285.9 

-133.4 

97.1 

2058 

6313 

-643.7 


Table 5.10b. Displacement, moment and membrane force data at points B and C. 
Solution by hierarchic shell model. Cartesian basis. Four finite elements. 


Pi 9 

N 

Me 

me 

me 

me 

me 

m B 

(M w )b 

(in) 

(lbf /in) 

(lbf/in) 

_0bf)_ 

_pw)_ 

(lbf/in) 

_(lbf)_ 

4,1 

224 

0.5363 

39276.0 

922.4 

150.3 

2108 

9420 

-704.5 

6,1 

456 

0.5415 

-2157.1 

-134.5 

94.0 

2055 

6318 

-645.4 

8,1 

784 

0.5418 

-283.8 

-133.2 

95.8 

2059 

6312 

-643.6 

10,1 

1208 

0.5420 

-285.5 

-132.7 

95.8 

2059 

6314 

-642.3 

12,1 

1728 

0.5422 

-283.7 

-132.5 

95.9 

2059 

6314 

-641.5 


- 36 - 



the cylindrical surfaces were exactly represented. In addition, there were slight 
differences in loading: In the three dimensional model half of the traction was ap- 
plied to the upper surface, half to the lower surface of the shell. In the hierarchic 
shell model the traction of 90 lbf/ft 2 was split between the upper and lower surfaces 
so that the force acting on the upper surface was equal to the force acting on the 
lower surface. 

Additional data are presented in Tables 5.10a,b. The definitions of N+, N v , 
A/*, M y are as follows: N+ and M x are defined on sections of constant <j>, see Fig. 



where <r+ is the normal stress in the r, <f>, y system; r< (resp. r 0 ) is the inner (resp. 

outer) radius of the shell and r m d = (r< + r 0 )/ 2 is the mean radius. Similarly N v and 
M v are defined on sections of constant y: 


Ny d = f <Ty dr, 

Jri 



(5.36) 


Convergence of the tabulated data is evident. When convergence is monotonic 
then some simple extrapolation scheme may be used to estimate limiting values. 
When convergence is oscillatory then the average of consecutive data is generally 
a good estimate of the limiting value. 


6. SUMMARY AND CONCLUSIONS. 


The only means for ensuring and verifying that engineering data computed 
from finite element solutions are within acceptable tolerance levels is by performing 
extensions. Extensions are orderly sequences of discretization constructed in such 
a way that the sequence of finite element solutions corresponding to a sequence 
of discretization converges to the exact solution. Of course, the exact solution 
and the norms in which convergence can be meaningfully measured depend on 
the choice of generalized formulation. In this paper the generalized formulation 
considered is the principle of virtual work. 

In conventional theories for curved beams, arches and shells the generalized 
formulation, as well as the displacement modes and stress distributions, vary from 
theory to theory. For this reason conventional theories do not constitute an orderly 
sequence, so that the corresponding solutions converge to the exact solution of a 
particular generalized formulation. 

In this paper we proposed an orderly sequence of discretizations controlled by 
two parameters. One of the parameters, denoted by p, represents the polynomial 
degree of basis functions with respect to standard coordinates mapped to the 
middle surface. The other parameter, denoted by q, represents the polynomial 
degree of the basis functions with respect to the standard coordinates which is 
mapped transversely to the middle surface. The finite element solutions converge 
to the exact solution of the fully three dimensional formulation of the theory of 
elasticity, based on the principle of virtual work, as p -» oo and q -+ oo. The 
examples presented herein are representative of a large class of problems which 
include thick, moderately thick and thin arches and shells. The results show that 
coarse meshes and g=l, 5<p<12, are generally sufficient to achieve levels of 
accuracy in the computed displacements, shear forces and moments which are 
normally expected in engineering practice. The hierarchic model characterized by 
q = 1 is similar to the Reissner-Mindlin theory except that one additional field 
which represents linear variation of the transverse displacement is incorporated. 
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The modeling strategy outlined herein views any beam, arch, plate or shell 
model as a particular discretization of the fully three dimensional model. Given a 
fixed finite element mesh, p-extension is performed when either p or q is increased. 
Convergence to the exact three dimensional solution is guaranteed. Recent surveys 
of the theoretical basis of p-extensions are available in [33,34]. The efficiency of 
p-extensions depends on the finite element element mesh. Initial mesh design 
is based on the analyst’s judgement concerning the expected smoothness of the 
solution. Modification of the initial mesh design is based on information generated 
by p-extension [24]. 

Our investigation has shown that curvilinear bases are not more advantageous 
than Cartesian bases in computer implementations. Originally, curvilinear bases 
had been used to permit analytical solutions for cylindrical, spherical, conical, 
toroidal and other shells mapped by relatively simple functions. In problems 
of current computer aided design and analysis a much wider range of mapping 
functions must be considered. Also, shells are connected to attachment lugs, 
stiffeners, etc. and are joined with other shells. In such regions the assumptions 
of conventional shell theories do not hold, and analytical treatment is impossible, 
yet those are the regions which are of greatest practical interest. The use of 
Cartesian bases provides for connection of shells of various types, stiffeners, three 
dimensional elements, etc., in a natural and systematic way. There is no need for 
transition elements and the process is computationally efficient. 

We have defined hierarchic sequences of basis functions based on Legendre 
polynomials. These basis functions lead to well conditioned stiffness matrices so 
that the accumulation of round-off error with respect to increasing p is slow. We 
have not encountered problems caused by round-off even at extreme aspect ratios 
(see Section 5.2) and p = 12. 

We have demonstrated by an example that in the case of thin arches and shells 
finite element discretizations can be substantially simplified by replacing E with 
E/( 1 - v 2 ) but otherwise using v = 0. The computed values of the stress resultants 
and displacements are not affected significantly by this substitution. 
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